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Abstract. We investigate the gravitational collapse of a spherically symmetric, 
perfect fluid with equation of state P = (T — l)p. We restrict attention to the 
ultrarelativistic ( "kinetic-energy-dominated" , "scale-free" ) limit where black hole 
formation is anticipated to turn on at infinitesimal black hole mass (Type II 
behavior). Critical solutions (those which sit at the threshold of black hole 
formation in parametrized families of collapse) are found by solving the system of 
ODEs which result from a self-similar ansatz, and by solving the full Einstein/fluid 
PDEs in spherical symmetry. These latter PDE solutions ("simulations") extend 
the pioneering work of Evans and Coleman (r = 4/3) and verify that the 
continuously self-similar solutions previously found by Maison and Hara et al for 
1.05 < r < 1.89 are (locally) unique critical solutions. In addition, we find strong 
evidence that globally regular critical solutions do exist for 1.89 < T < 2, that the 
sonic point for ~ 1.8896244 is a degenerate node, and that the sonic points 
for r > r dn are nodal points, rather than focal points as previously reported. We 
also find a critical solution for T = 2, and present evidence that it is continuously 
self-similar and Type II. Mass-scaling exponents for all of the critical solutions 
are calculated by evolving near-critical initial data, with results which confirm 
and extend previous calculations based on linear perturbation theory. Finally, we 
comment on critical solutions generated with an ideal-gas equation of state. 
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1. Introduction 

The formation of black holes is an exciting topic in general relativity, and a class 
of solutions which exists precisely at the threshold of black hole formation has 
recently generated considerable interest. These solutions have surprising properties, 
reminiscent of some thermodynamic systems near phase transitions, and, by analogy, 
have been called critical solutions. Critical phenomena in gravitational collapse were 
first discovered empirically in simulations of the massless Klein-Gordon field minimally 
coupled to gravity (EMKG) Q. Subsequent studies have shown that critical behavior 
is present in a variety of physical systems, and indicate that the phenomena are 
generic features of gravitational collapse in general relativity. In this paper we focus 
on the critical solutions for a spherically symmetric perfect fluid with equation of state 
P = (r — l)p, where p is the total energy density and T is constant, and present new 
solutions for T > 1.89. While we provide a brief introduction to critical solutions, 
the review by Gundlach is an excellent introduction to critical phenomena j2|, and 
additional information can be found in 101 . 
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1.1. Basic properties of critical solutions 

Imagine an experiment to investigate the details of gravitational collapse and black 
hole formation by imploding shells of fluid with a fixed equation of state. The initial 
energy density in the shell might be 



where Aq, ro, and A are parameters. In the course of the experiment we fix two of the 
three parameters, and allow only one of them, which we label p, to vary. For small p 
(assuming the fluid's initial kinetic energy is sufficiently large — i.e. in what one might 
call the ultrarclativistic limit), the fluid implodes through the origin and completely 
disperses. However, for p sufficiently large, in particular for p larger than some critical 
value p* , a black hole forms during the implosion, trapping some of the matter/energy 
within a finite radius. In the exactly critical case, p = p*, which represents the 
threshold of black hole formation, the evolution temporarily asymptotes to a special 
critical solution, Z*, which has a number of interesting properties, including scale 
invariance (self-similarity) and universality. The critical solution is universal in the 
sense that if we now use different "interpolating families" to probe the threshold of 
black hole formation, we will generically find the same critical solution (provided 
we remain in the ultrarclativistic regime). Additionally, in the super-critical regime 
p > p* , the black hole masses are well described by a scaling law 



Here the mass-scaling exponent 7 is also universal in the sense that it is independent of 
the particular choice of initial data family. (However, as first predicted by Maison 0, 
and Hara, Koike and Adachi J5], and as discussed in detail below, 7 is a function 
of the adiabatic index T.) 

One of the most profound consequences of the self-similar nature of critical 
collapse is that black hole formation in the ultrarelativistic limit turns on at 
infinitesimal mass. In analogy with second-order phase transitions in statistical 
mechanics, we refer to this behavior as Type II. As we will discuss shortly, Type I 
behavior, wherein black hole formation turns on with finite mass in interpolating 
families, has also been seen in various models of collapse and it is undoubtedly present 
in at least some of the perfect fluid models considered here. 

1.2. Critical solutions and one-mode instability 

A crucial feature of the critical solutions sketched above is that they are, by 
construction, unstable. If this is not obvious, one should observe that the critical 
solution is not a long-time (t — ► 00) solution of the equations of motion. Indeed, 
as sketched above, the only long-time stable "states" one finds from evolutions of a 
generic ultrarelativistic family of initial data either have (i) all of the fluid dissipated 
to arbitrarily large radii, with (essentially) flat spacetime in the interior, or (ii) some 
fluid dissipated to arbitrarily large radii, with a black-hole in the interior. The critical 
solution, Z* , on the other hand, exists just at the threshold of black hole formation, 
and, in near-critical evolutions, the dynamics asymptotes to Z* only during the strong- 
field dynamical epoch. For any given initial data, this strong-field regime persists 
for a finite amount of time (as measured, for example, by an observer at infinity). 
Eventually (and in fact, on a dynamical time scale) any non-critical data will evolve 
into one of the two stable end states. 
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Although the unstable nature of critical solutions was clear from the earliest 
phenomenological studies, considerable insight has been gained from the observation 
by Koike, Hara and Adachi |jj that the "sharpness" of the critical behavior seen 
in Type II collapse suggests that the critical solutions have exactly one unstable 
mode in perturbation theory. This ansatz immediately explains the universality of 
the critical solution: as p — > p*, one is effectively directly tuning out the single 
unstable mode from the initial data. Furthermore, using the self-similarity of the 
dynamics in the near-critical regime and a little dimensional analysis, it is an easy 
matter to relate the mass-scaling exponent to the Lyapunov exponent associated with 
the single mode. In fact, since the pioneering work by Koike et al , this picture of 
Type II critical solutions as one-mode unstable, self-similar "intermediate attractors" 
has been validated for essentially every spherically-symmetric model where Type II 
behavior has been observed in the solution of the full equations of motion. 

Moreover, the perturbative analysis applies equally well to Type I critical 
solutions which, arguably, have been well known to relativists and astrophysicists for 
decades, although perhaps not in the context of interpolating families. In this case, the 
critical solution is an unstable static or periodic configuration which, depending on how 
it is perturbed, will either completely disperse, or collapse to a finite-mass black hole. 
Once again, one generically finds that such solutions have a single unstable eigenmode, 
whose Lyapunov exponent is now a measure of the increase in lifetime of the unstable 
configuration as one tunes p — > p* . Type I behavior has been observed in the collapse 
of Yang-Mills [|| and massive scalar fields ||, and, as mentioned above, there is every 
reason to expect that it will occur in perfect fluid collapse. Indeed, it is well known 
that relativistic stellar models often exhibit a turn-over in total mass as a function of 
central density. Stars past this turn-over (generalized Chandrasekhar mass) are well 
known to be unstable, and, in fact, are almost certainly Type I solutions. 

1.3. Critical solutions and self- similarity 

Heuristically, systems exhibiting self-similarity appear identical over many different 
spatial and/or temporal scales, and generally arise in physical situations in which there 
are no natural length scales. Here it is important to note that a scale-free solution can 
be generated from a model which does have specific length scales, provided that the 
scaling solution represents a "self-consistent" limit. The current case of fluid collapse 
provides a perfect example. The rest mass of the fluid does set a length-scale, but 
the Type II critical solutions describe an ultrarelativistic limit wherein the rest-mass 
of the solution is irrelevant. To put this another way, we can have solutions of the 
equations of motion which have greater symmetry (scale symmetry in this case) than 
the equations of motion themselves. 

Self-similarity can be either continuous (CSS) or discrete (DSS), and both types 
have been observed in critical gravitational collapse. The perfect fluid critical solutions 
have continuous self-similarity of the first kind, a particularly simple self-similarity 
wherein the solutions can be written solely as functions of dimensionless variables, 
such as £ = r/t, where r is the radial coordinate in spherical symmetry and t is the 
coordinate time. An example of a CSS function of the first kind is shown schematically 
on a spacetime diagram in figure |]. 

At this point we should also note that the self-similar nature of Type II critical 
solutions provides a link between work on black-hole critical phenomena and the 
large body of literature dealing with self-similarity in the context of Einstein gravity 
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Figure 1. A schematic diagram showing a continuously self-similar (CSS) pulse 
at five different times as it moves toward the origin r = 0. The dotted lines are 
lines of constant f = r/t, the similarity variable. These lines converge at the 
space-time origin (r, t) = (0,0) in the upper left-hand corner of the plot, and the 
inset shows the pulse as a function of £. As the pulse moves toward the origin, it 
appears the same on smaller and smaller length scales. 

(see Carr and Coley |l(]] for a recent, extensive review). The self-similar ansatz 
has been widely employed, not only to produce more tractable problems, but also 
in investigations of possible mechanisms to generate counter-examples to the cosmic 
censorship conjecture. However, it is clear that not all self-similar solutions will be 
relevant to critical collapse, particularly if we restrict the definition of critical to "one- 
mode unstable" . Moreover, because most studies which are based on the self-similar 
ansatz have only considered the solutions themselves, and not perturbations thereof, it 
has proven non-trivial to identify which self-similar perfect-fluid solutions previously 
discussed in the literature are relevant to critical collapse Jb 



1.4- Review of previous perfect-fluid studies 

Shortly after the discovery of Type II behavior in scalar field collapse Evans, 
having been frustrated in attempts to "analytically" understand the massless-scalar 
critical solutions, turned his attention to perfect-fluid collapse. Armed with the 
intuition that the self-similarity of critical collapse was a defining characteristic, and 
aware of the existence of continuously self-similar relativistic fluid flows, he and 
Coleman jD| considered collapse in the context of the specific equation of state, 
P = gp. Significantly, they were able to construct a single critical solution, both from 
the self-similar ansatz (i.e. by solving ODEs), and by tuning the initial data used in 
solution of the full partial differential equations of motion. Moreover, they noted that 
a perturbation analysis about the inherently unstable critical solution would provide 
an accurate description of the near-critical dynamics, including the calculation of the 
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mass-scaling exponent 7. 

Such a perturbation analysis was quickly carried out (again for the case P — |p) 
by Koike et al Q], who, as mentioned above, made the crucial additional observation 
that the "sharp" transition in the mass scaling suggested that there was only one 
growing unstable mode associated with the critical solution, and that the Lyapunov 
exponent of the mode was simply the reciprocal of the mass-scaling exponent 7. Their 
analysis fully validated this conjecture — in particular, they found strong evidence 
for a single unstable mode with Lyapunov exponent 2.81055 . . . , corresponding to 
7 = 0.3558019 . . . , in excellent agreement with Coleman and Evans' "measured" value, 
7 re 0.36. 

At about the same time, Maison Q — assuming that the critical solutions would 
be continuously self-similar for other values of Y — adopted the CSS ansatz for the more 
general equation of state P = (T — l)p. He was able to construct CSS solutions for 
1.01 < T < 1.888, and additionally presented strong evidence that all of the solutions 
were one-mode unstable. Furthermore, the Lyapunov exponents, and hence the mass- 
scaling exponents were found to be T-dependent, with, 7, for example, varying from 
7 = 0.1143 for T = 1.01 to 7 = 0.8157 for T = 1.888. These calculations were 
particularly notable for providing early evidence that 7 was not a "truly" universal 
exponent, in the sense of having the same value across all possible models of collapse. 

One interesting outcome of Maison's linearized analysis of the equations of motion 
about the sonic point, was that at L ~ 1.888, two of the eigenvalues of the linearized 
problem degenerated, and the sonic point apparently changed from a node to a 
focus. This led Maison to conclude that regular self-similar solutions did not exist 
for r > 1.89. A similar analysis by Hara, Koike and Adachi || ^| (expanding on 
their previous work) again suggested a change in solution behavior at T sa 1.89. 
Those authors computed the CSS solutions, unstable eigenmodes, and eigenvalues for 
1 < r < 1.889, with results essentially identical to Maison's. Evans and Perkins [Q 
also performed the linear stability analysis for L < 1.888, finding the same results 
reported by Maison. In addition, they performed the first critical solution searches 
using the full set of PDEs for 1.05 < T < 1.5, confirming that the CSS solutions are the 
unique critical solutions for T is this range. Goliath et al 11 discussed, in the wider 
context of timelike self-similar fluid solutions, the mode structure of the P = (T — l)p 
CSS solutions, and reported that physical solutions do not exist for T > 1.89. More 
recently, after this paper appeared in pre-print form, Carr et al have extended this 
work for T > 1.89 |||l§. 

Furthermore, the conclusion of these linear perturbation analyses — that regular 
critical solutions for T > 1.89 do not exist — has inspired various proposals regarding 
the nature of T > 1.89 critical solutions || 10 ; |l7|| . (Recall that as long as we can 



set up interpolating data, there will be a critical solution, virtually by definition). 
One proposal is that a loss of analyticity at the sonic point for T > 1.89 violates a 
condition required to find the ODE solutions. Other proposals have suggested that 
the solution might become Type I, discretely self-similar, or display a mixture of DSS 
and CSS behavior. Some of these conjectures were evidently based on the fact that, 
under certain conditions, a stiff (P = p) perfect fluid can be formally identified with 
the EMKG system |]| p|, §(J. For example, it has been conjectured j| [l0|, that 
at some point as V —* 2, the critical solution might begin to display the discrete 
self-similarity characteristic of the EMKG critical solution. Brady and Cai [[t7] have 
previously computed threshold solutions for T < 1.98 using the full fluid equations 
of motion, finding — in all cases examined — evidence that the critical solutions are 
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both CSS and Type II. Using a two-step Lax-Wendroff numerical scheme to integrate 
the fluid equations, they calculated mass-scaling exponents by evolving supercritical 
initial data. However their code has severe resolution limitations, being able to observe 
scaling only over two orders of magnitude in \p — p* \ . 

Yet, lacking solutions for r = 2, and high resolution solutions near Y ~ 1.89, 
it was still expected that the perfect fluid critical solution changed its character as 
r — > 2. As we will discuss below, this does not seem to be the case, and in fact, 
r ~ 1.89 seems problematic only in the context of the the precisely self-similar ansatz. 
Specifically, we have strong evidence that the CSS ansatz generates an increasingly 
ill-conditioned problem as Y — > 1.8896244..., but that the PDEs remain perfectly 
well-behaved there. 

1.5. Summary of results 

Using a new perfect- fluid collapse code described in detail in pl| , this paper examines 
the perfect fluid critical solutions for the full range of Y, 1.05 < Y < 2, concentrating 
on those solutions for Y > 1.89. In addition to confirming the expected picture for 
1.05 < T < 1.889, we present further evidence that one-mode unstable CSS solutions 
exist in the regime Y > 1.89, up to and including, the limiting case r = 2. As part 
of this study, we are lead to re-investigate the CSS ansatz, and, in particular, the 
apparent change in solution behavior for Y ~ 1.889. In contrast to previous work, we 
solve the ODEs using arbitrary precision floating point arithmetic, which proves crucial 
to the analysis leading to our claim that the sonic points for Y > Y^n — 1.8896244 are 
nodal points, rather than focal points. We show that the CSS solutions for Y > 1.89 
are qualitatively identical to those for Y < 1.89, and verify that they can be computed 
from the full equations of motion (again using the code described in |2l| ). We also 
use simulations to compute mass scaling exponents in the regime Y > 1.89, where, 
interestingly, we find evidence to suggest that 7 — > 1 as Y — > 2. 

Finally, researchers in this field have tacitly assumed that it suffices to consider 
P = (r — l)p in the context of critical collapse, since (1) the critical solution 
naturally "drives itself" to the ultrarelativistic (kinetic-energy-dominated) regime and 
(2) P = (r — l)p is the only equation of state compatible with self-similarity |^2|, |23| . 
We explicitly verify this assumption (albeit on a case-by-case basis) by performing 
computations — using an ideal-gas equation of state — which display the anticipated 
behavior. 

2. Construction of the precisely critical solutions 

Spherically symmetric perfect fluid critical solutions can be constructed by searching 
for globally regular, spherically symmetric, CSS solutions to Einstein's equations p3| . 
As discussed in the Introduction, many studies of CSS perfect fluid spacetimes have 
appeared in the literature, and many properties of these solutions, such as their 
behavior near sonic points, are well known jl|, %, |o[ O, [2^, |4| ^5], ||(| ^7). The 
self-similarity ansatz reduces the equations to first order, autonomous ODEs, and 
this section focuses on the solution of these ODEs for Y > 1.89. We use the equations 
derived by Hara et al j| , and refer the reader to that paper for additional information 
concerning the derivation of the equations, and basic methods for their solution. 
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Figure 2. The geometric variable a of the critical solution for several values of 
r. The sonic point is at x = 0. The ODEs were integrated using the Maple V 
implementation of LSQDE with 30 digits and an absolute error tolerance e = 10 — 18 . 
We are increasingly unable to integrate these solutions outwards as T — > 2. This 
often occurs owing to a loss of numerical precision as u> — * 0, and the Lorentz 
factor, W = _ v 2 , becomes large (see figure [ij). 



2.1. The self- similarity ansatz 

We write the spherically symmetric line element in polar-areal coordinates as 

ds 2 = -a(r, if dt 2 + a(r, t) 2 dr 2 + r 2 (d9 2 + sin 2 9 ch/> 2 ) , (3) 

where the radial coordinate, r, directly measures the proper surface area. The Einstein 
equations give three relations for the metric quantities a and a (1: 16)— (1: 18) (equation 
numbers with a leading 'I,' refer to equations in pl[|), and there are two fluid equations 
of motion (1:23). A continuously self-similar spacetime is generated by a homothetic 
Killing vector £ |2£| , 

£i9ab = ^9ab- (4) 

We impose the homothetic condition and introduce coordinates s and x adapted to 
this symmetry 

s = -ln(-t), x = ln{~). (5) 

The time coordinate t is chosen such that the self-similar solution reaches the origin 
at t = 0, and the sonic point is at x = 0. We also define the dimensionless quantities 

N=— , A = a 2 , cj = 4vrr 2 a 2 p. (6) 
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Figure 4. The fluid velocity v of the critical solution for several values of I\ 
The sonic point is at x = 0. 
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Figure 5. The logarithm of the Lorentz factor for the fluid velocity, W = 
1/Vl — v' 2 , of the critical solution for several values of V. Notice the exponential 
growth in the Lorentz factor for larger values of I\ The sonic point is at x = 0. 



The self-similarity ansatz requires that all derivatives with respect to s vanish, thus 
reducing the system to autonomo us first order ODEs These equations are 

presented explicitly in [Appendix A ; here, we formally write them as 



M(y)y' = f(y), (7) 

where y is a "state vector" containing the four dependent variables y T — (A, N,v,u>), 
Ai is a 4 x 4 coefficient matrix, / is a "vector", and a prime (') denotes differentiation 
with respect to x. By construction, all solutions of the ODEs (Q) will be continuously 
self-similar; to find the unique critical solutions, we seek CSS solutions that are regular 
at both the origin and the sonic point |T^| . 

2.2. The sonic point 

The system of ODEs (Q) can be solved provided that the inverse matrix M. -1 exists, 
that is, provided that the determinant Q 

detMc< (1 + ^ )2 - (r - 1)(iV+ " )2 . (8) 
1 — v z 

is non-zero. When det M. = — a condition occurring at a sonic Cauchy horizon 
or sonic point — the ODEs cannot be integrated without further assumptions. In 
particular, if detA^ = 0, the derivatives, y' , may either (i) not exist or (ii) be 
undefined. In the former case, the functions may be continuous but not differentiable, 
or a shock may form (discontinuous functions) at the sonic point; the latter case 
corresponds to the physically-relevant regular solutions in which we are interested. 
By definition, the sonic point is the position where the magnitude of the fluid velocity, 
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as measured by an observer at constant x — ln(— r/T), is equal to the fluid sound 
speed 

c s = Vr^l. (9) 

All CSS perfect fluid solutions that are regular at the origin have at least one 
sonic point |22j. If y is regular at the sonic point, then the rows of M. must be linearly 
dependent so that det M. = 0. This condition allows one to parameterize the CSS 
solutions with a single parameter, t> sp , which is the fluid velocity at the sonic point. 
Furthermore, this regularity condition fixes y' sp in terms of v sp and v' sp , yielding a 
quadratic condition on v' sp which we schematically write as 

c 2 ^p 2 + civ' sp + c = 0. (10) 

Here, the coefficients, cq, c\ and C2, are complicated functions of y ap , and for simplicity 
of presentation will not be given explicitly. The key point is that this constraint — 
that the critical solutions are regular at the sonic point — limits the number of solutions 
to discrete values of v sp , and virtually eliminates the possibility that globally regular 
solutions with more than one sonic point exist ^7). Indeed, all of the T < 2 critical 
solutions we have found, either from a CSS ansatz, or by solving the full Einstein/fluid 
equations, have only one sonic point. 



2.3. Solving the ODEs 

The system of ODEs (|7|) is solved by choosing a candidate fluid velocity, v sp , at the 
sonic point, and integrating numerically from the sonic point toward the origin. The 
inward integration is halted when either A < 1 or detAi = 0, and these generic 
stopping criteria allow one to determine the parameter v sp by a standard "shooting" 
procedure. (If A < 1 signals that v sp is too small, then det M = indicates that v sp is 
too large, and vice versa). Once v sp has been determined from the inward integration, 
the solution can be completed by integrating outwards from the sonic point. The entire 
solution process is complicated by the fact that the integration can not actually begin 
at the sonic point, since det.M = there. Therefore, we first expand the dependent 
variables y about the sonic point to first order 

Vo ~ y sp + y' sp Ax, (11) 

where Ax = x — x sp , and actually begin the integration from x Q . Ax is chosen so 
that the 0((Ax) 2 ) error terms in the expansion, are smaller than the error tolerance 
allowed in the solution. We obtain y' sp by solving (jl^) for v' sp and integrate the ODEs 
for both roots. 

The ODEs are integrated using LS0DE, a robust numerical routine for integrating 
ODEs [^9], [50), and all of the critical solutions can be found using double precision 
arithmetic, except those solutions for T sa 1.89. These T 1.89 solutions require 
greater precision, and can be found using the arbitrary precision implementation of 
LS0DE in Maple V, which also proved invaluable for convergence testing the solutions. 
In the convergence tests we vary Ax, the LS0DE absolute error tolerance e (the relative 
error tolerance is set to zero), and the number of digits used in the calculation, while 
monitoring the residual of the algebraic constraint (A. 6) as an indication of the error 
in the solution. For example, we calculated the critical solution for T — 1.99 using 40 
digits and error tolerances e = 10~ 10 , 10~ 15 , 10~ 20 and 10 -25 , and then performed 
similar tests using 30 and 35 digits. In all cases the solutions converge, and the residual 
of (O) is O(e). 
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The critical solutions for several values of T are shown in figures One notes 
that as r — > 2, it becomes increasingly difficult to integrate outwards from the sonic 
point, and this limitation is especially apparent in the Lorentz factors for the T = 1.95 
and r = 1.99 critical solutions shown in figure ||. The ODE solver fails when the 
Jacobian matrix, 

becomes ill-conditioned, which often occurs when w»£ and/or v f=a 1 — s, where e is 
the error tolerance supplied to LSODE. We can extend these solutions further outwards 
simply by increasing the number of digits used in the calculation. However the increase 
in maximum physical radius as a function of resolution is very modest — since the 
Lorentz factor W increases exponentially with increasing radius — and the process 
quickly becomes prohibitively expensive. Clearly, J is always poorly conditioned 
near sonic points, but J is particularly poorly conditioned near the sonic point for 
T sa 1.89 — one can view this as a clear signal that numerical work will be difficult in 
the r w 1.89 regime. 

2.4- Nature of the sonic point for T > 1.89 

Although nonlinear systems of ODEs are often impossible to solve in closed-form, 
qualitative features of their solutions can frequently be deduced by linearizing the 
equations about "critical" points. Perfect fluid CSS solutions have often been studied 
using this type of analysis |fl 24, |25|, 26, and here we discuss some of 



these results in the context of our work. We emphasize, however, that we have not 
performed perturbation analyses in our current work. 

The CSS perfect fluid equations (Q) can be linearized about the sonic point, 
resulting in a system we can write in the form 

y' = By, (13) 

where B is a matrix which, gcncrically, has two non-zero eigenvalues, which we label k\ 
and K2 (or simply k when the distinction is irrelevant), with corresponding eigenvectors 
Vi and V 2 - Near the sonic point, the solution of the linear equations (|l3|) can be 
written [^5| 

y = y sp + + k 2 V 2 e K2X , (14) 

where k± and k 2 are arbitrary constants. The eigenvalues K provide important 
information about the solutions near the sonic point, and we classify the sonic point 



according to the relative values of K, as given by the quantity |22 



*=^. (15) 



Here, 



and 



77 = 4(3r - 2)U 2 - (3r 2 - 12r + 8)(1 - AU), (16) 

"-^ (17) 

The sonic point classification in terms of $ is shown in table Due to the facts 
that (i) i9 is only a function of y sp , and (ii) the eigenvalues n are related to v' sp p7| , 
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Figure 6. The sonic point of the T^n — 1.8896244 critical solution is a 
degenerate node, and this figure shows the approach to degeneracy by plotting 
the quadratic function in (llCj) as V — > r c i n . The roots of QIOI) are possible values 
for v' at the sonic point, and when T = r<j n these roots are equal, i/, . Each 
frame shows a parabola for three values of T. The parabolas on the left are for 
r less than r<j n (r<), and the parabolas on the right are for T greater than 
(r>). The center parabola has the same T in all five frames, namely Vq — ^dn 
(Pc = 1.88962441796875), and one must be careful not to mistake it for a vertical 
axis in frames (a)— (c). For all T, the critical solution's v' sp is the root closest to 
the center (v' d J. We estimate T dn ~ 1.8896244169921874 (v Bp = -0.18696..., 
v' Bp = 1.7385 . . . , and i3 sp = 1.0000000002), and for clarity, we use v' dii to translate 
the horizontal axis such that the parabolas cluster around zero, and normalize the 
parabolas. In frame (a), T< = 1.889 and T> = 1.890. In (b), T< = 1.889624 
and T> = 1.889625. In (c), T< = 1.88962425 and T> = 1.8896245. In (d), 
T< = 1.889624375 and T> = 1.8896244375. In (e), T< = 1.8896244140625 and 
T> = 1.889624421875. These calculations were done with MapJe Vusing 30 digits 
and e = 10 -18 . 
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Table 1. Classifications of the sonic point using 

i? < saddle point 

i? > nodal point 

•d £ C focal point 

i9 = 1 degenerate nodal point 



we can make a connection to the linearized theory without explicit calculation of the 
eigenvalues. 

Maison and Goliath et al [|ll] have previously concluded that the sonic points 
for r > 1.89 are foci, with complex k and v' sp , and hence have suggested that physical 
self-similar solutions do not exist for T > 1.89. (Hara et al did not address the 
existence of solutions for T > 1.89, but presumably also encountered problems with 
their numerical analysis in that regime.) However, we find that i? > for all T > 1.89 
critical solutions, and thus conclude that the sonic-points for those solutions are 
nodes rather than foci. It seems plausible that this apparent contradiction stems 
from insufficient numerical precision in the earlier studies. To provide some specific 
evidence to back this claim, we have used Maple with 30 digits to find a critical 
solution for T ~ Tdn- Then, taking v sp from this solution, we have calculated •& using 
both 30 digits in Maple and FORTRAN double precision. The FORTRAN calculation gave 
a complex ■& — which would support the erroneous (we claim) conclusion that the sonic 
point is a focus. The same calculation done with greater precision using Maple shows 
that the sonic point is actually a node. In addition, we find that for T < T^n (with 
r restricted to T > 1.8 for simplicity, and T^n defined below), the critical solution's 
v' sp is the maximum root of (|To|), while for T > Tdn, the critical solution's v' sp is the 
minimum root. As T — > Tdn, the two roots v' sp come closer together until they are 
equal for Tdn, as shown in figure [| Here the sonic point is a degenerate node with 
r d „ ~ 1.8896244 { Vdn = 0(e), and tf d „ = 1 + O^/i)). 



3. Simulation of critical solutions 



A crucial check that the CSS solutions of the ODEs are indeed the unique critical 
solutions involves a comparison with the solutions of the full Einstein/fluid equations. 
Although relativistic fluid equations are known to be difficult to solve (particularly 
relative to "fundamental-field" equations, such as the EMKG system), modern high- 



resolution shock-capturing schemes [31, |32|, 33, |3J] allow one to calculate flows with 
shocks of almost arbitrary strength. However, perhaps the greatest challenge in finding 
the perfect fluid critical solutions — especially as T — > 2 — is the accurate treatment of 
flows with very large Lorentz factors. In [2l| we have outlined our computer program 
which solves the spherically-symmetric Einstein/fluid system, and using this code we 
have found the perfect fluid critical solutions for 1.05 < T < 2. All of these solutions 
are continuously self-similar (CSS) and black hole formation for near-critical initial 
data begins with infinitesimal mass (Type II). In this section we compare the ODE 
and PDE solutions, discuss the mass-scaling exponents 7, and finally, briefly discuss 
critical solutions for the ideal-gas equation of state. 
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Figure 7. The T = 1.9 critical solution. The solid lines are the solution obtained 
by solving the ODEs, and the triangles indicate selected points from the solution 
of the PDEs. 



3.1. Critical solutions 

Figures f?]-|| show the PDE critical solutions for T = 1.9 and T = 1.99 respectively, 
and include the ODE solutions for comparison. The r = 1.9 PDE and ODE solutions 
compare well, while figure || indicates that the T = 1.99 PDE solution underestimates 
a and uj. From our experience with other critical solution searches, we feel that this 
discrepancy is the result of insufficient (spatial) resolution. As T — > 2, the fluid 
becomes increasingly dynamic, requiring greater precision to resolve the solution's 
relevant features, and it becomes increasingly expensive to calculate the critical 
solutions. The mass-scaling exponents 7 shown in table |^ provide a quick guide to the 
requisite dynamical range for a critical evolution as T — > 2. As (p — p*) changes by 
n orders of magnitude, the relevant length scales in the solution, such as the radius 
of a black hole i?BH, change by "fn orders of magnitude. The mass-scaling exponent 
for the stiff fluid, 7 w 1, is almost three times larger than the scaling exponents for a 
radiation fluid (7 « 0.36) or massless scalar field (7sf ~ 0.37), and simulations of the 
stiff fluid critical solutions require correspondingly more resolution. 

The critical solutions for T < 2 all appear very similar; indeed one can imagine 
that one could smoothly transform a solution for a given Y into a solution for a different 
r. At first glance, the T = 2 solution seems to fit nicely into this "family" of critical 
solutions parameterized by T — it is CSS, Type II, and differs only slightly from the 
r = 1.99 critical solution. However, the ODE solution (obtained by solving ODEs with 



Critical phenomena in perfect fluids 



15 



_ 1 1 1 1 1 1 1 1 1 1 1 1 1 

: a / 


I i i i i i i _ 

Aa a a = 

A — 






z CO A 




*± = ' A A A 


V 


AAA A 

A - 
A — 

i 

— 


Z i 

A * ] 


r = 1.99 -_ 


i i i i i i i i 1 i ifSrt 


i i i i i i i 



-3-2-10 1 2 

x 



Figure 8. The V = 1.99 critical solution. The solid lines are the solution obtained 
by solving the ODEs, and the triangles are selected points from the PDE solution. 
The PDE solution underestimates the fluid density in the pulse leading, to a 
corresponding error in a. This problem stems from a lack of resolution in the 
computational grid. 



the CSS ansatz) indicates that important differences may exist between the r = 2 and 
r < 2 critical solutions. As noted previously, we are unable to integrate the ODEs 
for r near 2 to arbitrarily large x. In these cases, we observe that the Lorentz factor, 
W, grows exponentially (see figure [5]), with a corresponding exponential decrease in 
u>, until LSDDE is unable to satisfy the required error tolerances. We emphasize that 
in these T < 2 solutions, the fluid velocity and density retains its expected "physical" 
properties: u> > 0, and |t>| < 1. The r = 2 solution (see figure ||), on the other hand, 
displays very different behavior. Instead of the exponential approach of v — > 1 and 
uj — > 0, we find that u) and v pass through their expected physical bounds, giving 
us < and v > 1. As is generally known, the stiff perfect fluid can be related to a 
scalar field, and, very recently, a CSS scalar field solution has been found by Brady 
and Gundlach that exactly matches the ODE T = 2 fluid solution p5] . 

The r = 2 PDE solution is also shown in figure ^|, and in contrast to the ODE 
solution, this solution retains "physical" values for the fluid variables. However, this 
property of the PDE solution is achieved by fiat: we impose a "floor" on the fluid 
variables such that p > and |?;| < 1 plf - While the floor is used generally for 
r > 1.8 without noticeable ill effect for T < 2, it is clear that the floor affects the 
r = 2 solution (as compared to the ODE solution) even in the regime where u> > 
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Figure 9. The Y = 2 critical solution. The solid lines show the solution obtained 
from the ODEs, and the dotted lines with triangles show the solution obtained by 
solving the PDEs. Here some divergence between the PDE and ODE solutions 
can be seen. The PDE solution for uj lies above the ODE solution for x < 0, 
behavior opposite from that observed in the T = 1.99 solution. Beyond the 
sonic point, x = 0, the solutions become very different. Here the ODE solution 
becomes unphysical — in the traditional understanding of a perfect fluid — with 
lu < and v > 1. The PDE solution retains physical values for v and u) — v < 1 
and lo > — because these values are constrained to the physical regime in the 
evolution code 1211 . 



and \v\ < 1. The mathematical and physical significance of this observation is clearly 
an issue which requires more study — for example, can a vacuum region be matched to 
the r = 2 fluid in such a way that the fluid remains equivalent to a EMKG field |3(| ? 

Finally, we note that the T = 2 critical solution is not related to other familiar 
EMKG solutions, such as the Roberts solution ^8), or the EMKG critical 
solution M. Both of these solutions have spacelike gradients of the scalar field. 
Extracting data from a near-critical T = 2 solution, we have set equivalent initial 
data for an EMKG evolution and then have evolved the data with the Einstein/scalar 
equations of motion. However, the evolution of the scalar field does not seem to match 
the r = 2 perfect-fluid critical solution for any appreciable length of time, and naive 
variations of the initial data apparently produce the usual DSS scalar field critical 
solution at the black-hole threshold. 
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3.2. Mass-scaling exponents 

Mass-scaling exponents 7 are found by evolving near-critical initial data sets which 
lead to the formation of black holes. In our coordinate system, black hole formation 
is signaled by 

2m(r, t) 



-1, (18) 

where -Rbh is the (areal) radius of the black hole. The black hole mass is then simply 
given by 

Mbh = ^p- (19) 

As mentioned earlier, all of the critical solutions discussed here are Type II, meaning 
that the associated black-hole transition begins with infinitesimal mass. As a typical 
example of our results, the mass-scaling of near-critical solutions for T — 2 is shown 
in figure [lC]. 

The simple adaptive grid that we use Q did not allow us to calculate Mbh with 
sufficient accuracy to justify searching for p* to the limit of machine precision in a 
reasonable amount of time. We therefore estimated p* by searching for the best linear 
fit to 

In Mbh oc j\n\p~p*\ . (20) 

The totality of mass-scaling exponents 7 calculated from our simulation data are shown 
in table ||, along with the values predicted from Maison's perturbative calculations Q . 
(These exponents for Y < 1.98 are similar to those found independently by Brady and 
Cai ]l7j.) For a variety of reasons, estimation of the error (no doubt overwhelmingly 
"systematic" ) in the mass-scaling exponents is not an easy task, and we therefore have 
provided estimates of 7 which are conservative in their use of "significant" digits. One 
notes that as Y — > 2, the mass-scaling exponent 7 — > 1. This trend suggests that 
7 = 1 for r = 2. However our numerical results can only determine 7 to very limited 
precision, and we are not currently aware of any argument for a precise equality. 



3. 3. Critical solutions for the ideal gas 

The equation of state (EOS) P = (T — l)p can be interpreted as the ultrarelativistic 
limit of the ideal-gas EOS 

P=(Y-l)p e, (21) 

where p is the rest energy density and e is the specific internal energy density. 
Following Ori and Piran p2| , Evans (2^| (and others), we have argued |Q that self- 
similar perfect fluid solutions require the ultrarelativistic EOS. Let us now consider 
searching for critical solutions with the ideal-gas EOS. Heuristically, one can describe 
critical behavior in terms of competition between the fluid's kinetic energy and 
gravitational potential energy. One might expect that in the critical solution, which 
stands just on the verge of black-hole formation, P ^> p a , and that the critical solutions 
for the ideal-gas EOS would correspond to the ultrarelativistic EOS solutions. Using 
an evolution code, we have found critical solutions for the ideal-gas EOS, and these 
solutions do match the corresponding solutions with the ultrarelativistic EOS. As an 
example, the critical solution for a Y — 1.4 ideal gas is compared with the precisely CSS 
ultrarelativistic Y — 1.4 solution in figure [D]. Additional evidence that near-critical 
ideal gas solutions are ultrarelativistic is shown in figure O. 
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Figure 10. Illustration of black-hole mass scaling for the case T = 2. In this 
instance — as for all values of V considered here — the critical behavior is Type II, 
allowing one to create arbitrarily small black holes through sufficient fine-tuning 
of initial data. 



Table 2. The mass-scaling exponent 7 as a function of the adiabatic constant 
r. The second column shows the mass-scaling exponents estimated from A/bhCp) 
by evolving near-critical (p — » p*) initial data. For comparison, "calculated" 
exponents — computed from perturbative calculations — are also listed. 
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Figure 11. A comparison of the critical solutions for the ideal-gas and 
ultrarelativistic equations of state for T = 1.4. The solid lines are solutions 
obtained by solving the ODEs with the ultrarelativistic EOS, and the triangles 
are selected points from the PDE solution which uses the ideal-gas EOS. Th p 
ideal gas is in the ultrarelativistic limit near the infalling matter (see figure 120, 
and the two solutions correspond in this region. At large r the ultrarelativistic 
approximation breaks down, and the solutions differ. The ideal-gas EOS solution 
was computed using a code similar to the one described in [ElJ . 



4. Conclusion 

Following seminal work by Evans and Coleman jlj|, Maison and Hara, Koike 
and Adachi || have previously reported the existence of CSS solutions with a single 
unstable (relevant) mode for T < 1.89. In this paper we have shown that such 
solutions also exist for T > 1.89, and have found evidence that the sonic point is 
a degenerate node for the Tdn — 1.8896244 critical solution. Our results come from (i) 
a demonstration of the existence of globally regular CSS solutions for T < 2 (starting 
from a CSS ansatz), and (ii) the evolution of finely-tuned initial data in full simulations 
of the Einstein/fluid equations. The simulations verify that the CSS solutions sit at 
the threshold of black hole formation, and also allow us to compute mass scaling 
exponents which are in good agreement with predictions from perturbation theory. 
The r = 2 critical solution is also CSS and Type II, with a mass-scaling exponent 
very close to (if not identical to) unity. We have also investigated critical collapse 
using the ideal-gas equation of state, and, as expected, have found that the fluid is 
well-described in the near-critical regime by an ultrarelativistic approximation, and 
thus has critical solutions identical to those generated using P = (T — I) p. 

While we have addressed some of the questions regarding the T > 1.89 critical 
solutions that have appeared in the literature, other avenues for further research 
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Figure 12. P and p a are shown at the moment of maximum compression in this 
log-log plot of a marginally subcritical evolution of a Y = 1.4 ideal gas. Note that 
near the origin, the ideal gas is clearly in the ultrarelativistic limit where P>p„. 
At large r (X), the ultrarelativistic limit no longer holds, and the solution does 
not match one computed using the ultrarelativistic EOS. 

remain. Perhaps foremost is the need for a greater understanding of the T = 2 
critical solution, both in a more precise estimate of its mass-scaling exponent, and its 
precise relation to the EMKG field. In addition, Gundlach |3!| and Martin-Garcia 
and Gundlach jl0| have recently investigated the critical behavior of an axisymmetric 
radiation fluid with angular momentum. Full simulations of such scenarios will 
undoubtedly yield further interesting results. 

Finally, we wish to point out an obvious use of critical solutions for testing and 
benchmarking new codes which solve the equations of motion (PDEs) for general 
relativistic fluids. In particular, we are unaware of widely-used test-bed solutions 
which combine both dynamic fluid flows and strong gravitational fields. Common 
tests for such codes currently include Riemann shock tube configurations, static 
Tolman-Oppenheimer-Volkoff solutions and Oppenheimer-Synder dust collapse. The 
former tests lack gravitational effects, and the latter tests lack either dynamic fluid 
flows or pressure effects. While we do not advocate abandoning traditional tests, 
we wish to emphasize that critical solutions are novel in that they combine several 
computationally challenging characteristics: extremely relativistic and dynamic fluid 
flows, strong gravitational field dynamics, and relevant features over many space/time 
scales. In particular, as adaptive mesh refinement (AMR) becomes available for multi- 
dimensional simulations, these features should make critical solutions ideal candidates 
for code calibration. 
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Appendix A. ODEs for a self-similar spacetime 

The Einstein equations for a spherically symmetric, CSS perfect fluid are presented in 
this appendix for reference. These equations are derived in ||. 

A' , 2w(l + (r-l)i; 2 ) lk x 

A =1 - A+ 1-,* ^ 

^ = -2 + A-(2-T)oj (A.2) 



A' ZTNvlo 



A 1 



(A.3) 



lo 1 — ir 

= -(2-r)iVi;- ^—ANv + (2-T)Nvlo (A.4) 

.u)' r(l + Nv)v' 

7V — fi 9 — 3F 
= (2-T)(T-l)Nu + ^— -N + ^^-AN (A.5) 



The equations ( A.2 )-( |A.5[) are integrate d to find th e cr itical solutions, and comprise 
the system My' = f in (|7j). Equations ( |A.l| ) and (A..3) can be combined to give an 
algebraic constraint 

(1 - A)(l - v 2 ) + 2lj (1 + (r - l)v 2 ) + 2TNvlj = 0. (A.6) 

This equation is not used to integrate the ODEs, but is monitored during the 
integration as an indication of the error in the solution. 



References 



[l] 

[2] 
[3] 

[5] 
[6] 
[V] 



Choptuik M W 1993 Phys. Rev. Lett. 70 9 

Gundlach C 1998 Ad v. Theor. Math . Phys. 2 1; ;r-qc/9712084 



gr-qc/9803075 



Choptuik M W 1998 

Maison D 1996 Phys. Lett. B 366 82 

Hara T, Koike T, and Adachi S 1996 gr-qc/9607010 

Koike T, Hara T, and Adachi S 1999 Phys. Rev. D 59 

Koike T, Hara T and Adachi S 1995 Phys. Rev. Lett. 74 5170 



Critical phenomena in perfect fluids 



22 



Choptuik M 
Brady P 1997 




Bizon P 1996 Phys. Rev. Lett. 77 424 



Carr B J and Coley A A 1998 gr-qc/9806048 

Goliath M, Nilsson U S and U ggla C 1998 Cla ss. Quantum Grav. 15 2841 
Carr B J and Coley A A 1999 |gr-qc/9901050 
Evans C R and Coleman J S 1994 Phys. Rev. Lett. 



72 1782 

Perkins T J W 1996 thesis (unpublished) 

Carr B J, Coley A A, Goliath M, Nilsson U S and Uggla C 1999 gr-qc/990131 
Carr B J, Coley A A, Goliath M, Nilsson U S and Uggla C 1999 gr-qc/990270 
Brady P R and Cai M J 1998 ;r-qc/9812071 
Taub A H 1959 Arch. Rat. Mech. Anal. 3 312 
Madsen M K 1985 Astrophys. Space Sci. 113 205 
Madsen M K 1988 Class. Quantum Grav. 5 627 

Neilsen D W and Choptuik M W 1998 submitted to Class. Quantum Grav. 
Ori A and Piran T 1990 Phys. Rev. D 42 1068 

Evans C R 1993 in Lecture Notes of the Numerical Relativity Conference, Penn State University, 

1993 (unpublished) 
Bogoyavlenskii O I 1977 Sov. Phys. JETP 46 633 
Bicknell G V and Henriksen R N 1978 ApJ 219 1043 
Bicknell G V and Henriksen R N 1978 ApJ 225 237 
Foglizzo T and Henriksen R N 1993 Phys. Rev. D 48 4645 
Cahill M E and Taub A H 1971 Comm. Math. Phys. 21 1 

Hindmarsh A C 1983, in Scientific Computing Stepleman R S et al (eds.) (North-Holland, 

Amsterdam) 55 
Petzold L R 1983 J. Sci. Stat. Comput. 4 136 

LeVeque R J 1992 Numerical Methods for Conservation Laws (Birkhaiiser-Verlag, Basel) 
LeVeque R J 1998 , in Computational Methods for Astrophysical Fluid Flow, 27th Saas-Fee 
Advanced Course Lecture Notes (Springer- V erlag. Berlin, to be published); also available at 



ittp: //sirrah. astro .unibas . ch/saas-f ee/ 



Ibanez J M— , Marti' J M— , Miralles J A and Romero J V 1992 Approaches to Numerical Relativity, 

(Cambridge University Press, Cambridge) p 223. 
Romero J V, Ibanez J M^, Marti' J M& and Miralles J A 1996 ApJ 462 839 
Brady P, Gundlach C, Neilsen D, and Choptuik M 1999 in preparation 
Gundlach C 1999 personal communication 
Roberts M D 1985 Gen Rp.L Cray 17 913 
Roberts M D 1998 



-qc/981109: 



Gundlach C 1998 Phys. Rev. D 57 7080 
Martin-Garcia J M and Gundlach C 1998 |gr-qc/9809059 



